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ABSTRACT 

The  paoer  solves  the  HMMS  disaggregation 
model  using  numerical  methods.  The  first  oroDOsed 
soproach  is  to  turn  the  original  problem  into  en 
unconstrained  minimization  and  apoly  lattice- 
Fibonacci  search.  Then  the  model's  solution  is 
considered  as  a  root-finding  problem  and  two 
approaches  are  compared:  Bolzano's  method,  and  a 
modified  lattice-Fibonacci  technique  proposed  in 
the  paper.  These  search  algorithms  proved 
extremely  efficient  in  yielding  solution  estima- 
tes with  errors  of  the  order  of  tenths  of  a  per- 
cent and  less  in  a  small  number  of  iterations. 
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1.  INTRODUCTION 

The  issue  of  making  aggregate  decisions  /I/  in  caoacity 
olanning  nroblems  is  a  direct  consequence  of  the  uncertain  environraeni 
in  which  oroduction  takes  piece  and  of  our  limited  ability  in 
gathering  and  processing  very  large  amounts  of  detailed  data. 
Moreover,  a  manager's  approach  to  the  capacity  planning  question 
is  by  its  nature  aggregate  in  the  first  place,  in  order  to  give 
him  a  broad  view  of  where  he  stands.  Then,  naturally,  the  issue  of 
disaggregating  the  information  generated  by  the  aggregate  analysis 
stage  is  the  next  question,  which  is  in  no  way  simpler  than  the 
first  one  /2/. 

The  book  of  Holt,  Modigliani,  Muth  and  Simon  /3/  is  a 
self-contained  work  that  addresses  both  decision  processes  pointed 
out  above.  By  using  a  quadratic  inventory  cost  function  together 
with  the  overtime,  idle  time,  hiring,  layoff,  and  other  related 
cost  functions  (all  quadratic),  the  authors  obtain  linear  decision 
rules  for  making  optimal  aggregate  decisions  for  each  period. 
The  aggregation  is  an  extreme  one  of  all  product  types  into  a 
single  category  requiring,  of  course,  the  use  of  appropriate  com- 
patible units  that  allow  the  transformation  to  be  made.  The  linear 
decision  rules  specify  the  aggregate  oroduction  and  work  force  for 
each  period.  Then,  using  the  aggregate  decision  as  a  constraint  on 
the  other  more  detailed  and  numerous  decisions  that  concern  the 
production  of  individual  items,  rules  are  obtained  to  yield  the 
optimal  decisions  under  the  constraint. 

There  are  4  disMggregation  models  presented  in  the 
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HMMS  book: 

1-  determining  oroduction  quantities  by  minimizing 
inventory  holding  costs  and  set-up  costs  subject  to 
the  aggregate  inventory  constraint  (/3/,ch.  10); 

2-  determining  buffer  stocks  by  trading  off  inventory 
costs  for  stock-out  costs  under  the  aggregate  con- 
straint end  disregarding  the  batch  characteristic  of 
the  production  (/3/,  ch.ll); 

3-  planning  order  points  given  fixed  lot  sizes,  by 
minimizing  inventory  holding  costs  and  inventory 
depletion  costs  under  aggregate  constraint  (/3/,  ch.l2); 

4-  determining  lot  sizes  and  safety  stocks  by  minimizing 
setup  costs,  inventory  holding  and  depletion  costs, 
under  aggregate  constraint  (/3/,  ch.l3). 

Models  2,  3,  and  4  above  require  the  estimation  of  the 
cost  of  being  out  of  stock  by  one  unit  in  order  to  comoute  the 
inventory  depletion  costs.  From  the  managerial  ooint  of  view  this  is 
a  difficult  task,  hard  to  implement  because  the  intangibles  invol- 
ved in  stock-outs  can  not  always  be  cast  in  precise  mathematical 
forms.  This  is  why  the  approach  of  computing  safety  stocks  and  order 
points  based  on  service  level  considerarions  /4/ch.6,  /5/  seems 
to  be  more  practical  and,  possibly,  with  more  managerial  apoeal. 

This  paper  proposes  to  analyze  model  1  above,  from  the 
computational  viewpoint,  to  find  an  efficient  way  to  solve  it 
using  numerical  methods,  and  then  compare  the  results  with  the 
solution  given  by  the  authors. 


where: 
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2.  HMMS  MODEL  FOR  COMPUTING  LOT  SIZES 
UNDER  AGGREGATE  CONSTRAINT 

2.1.  The  model 


s.t. 

n_   Q^ 


(1) 


i   a  denotes  item  i,  i=l,2,...,n; 

n 
Fi  =  setuD  coat  for  a  lot  ; 

Si   «  forecasted  sales  rate  in  units  oer  period  of  time; 

Q,   =  lot  size  in  units  of  product; 

Cji    =  cost  of  holding  one  unit  of  inventory  one  period 

of  time; 
n,   =  factor  for  converting  units  of  the  i'th  oroduct  to 

the  corresponding  number  of  common  units; 
I^  =  aggregate  inventory  in  the  common  unit. 

Since  the  purpose  of  the  paper  is  to  study  the  solution 
to  the  model,  it  was  left  in  the  form  oreeented  by  the  authors. 
It  is  clear  that  if  we  also  consider  the  problem  of  safety  stocks 
the  constraint  in  (1)  should  be  expressed  in  terms  of  aggregate 
production  P: 

S  u^Qi  =  P  (2) 

rather  than  Iq.  This  is  because  aggregate  cycle  stock  cannot  be 
computed  unless  we  know  t)ie  aggregate  safety  stock,  which  in  turn 
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can  be  computed  (under  service  level  considerations)  only  after 
Qj  8  are  known.  Fortunately  the  method  of  solving  (1)  remains 
unchanged  if  (2)  becomes  the  constraint. 

By  using  the  Lagrangian  multiolier  method  in  (1)  we 
tranaforme  the  constrained  minimization  into  an  unconstrained  one, 

min  T  -^  r  ^Fi^i  ^  CjiQi  n   u.Q, 

min  L  =>_^  (  Q    ■»■  — 2 —  '   "*"^^^Q~  ^^   "5 — ^ 

The  first-order  conditions  for  a  minimum  are  ~^  =  0 

yielding: 

■2CpiSi 


^li  ""i 

Replacing  Qj  into  the  constraint  in  (1)  leads  to  an 
equation  whose  only  unknown  is  A  : 


I  2C„.S.' 


Iq  (4) 


Since  n  is  usually  large,  equation  (4)  is  imoossible  to 
be  solved  exactly  for  A . 

2.2.  Authors'  solution 

We  are  interested  here  in  the  general  case  of  (3),  when 
no  special  simplifying  assumotions  are  made  such  as  :  constant 
sales  composition,  identical  holding  costs. 

The  first  solution  orooosed  by  the  authors  is  a  gra- 
phical method.  By  drawing  a  smooth  eye-fitted  curve  through  a  few 
plotted  Dointa,  one  can  obtain  the  grapij  of  the  relation  between 
Iq  and  A  for  any  fixed  set  of  salesrates  S,,S2f...,S  ,  Given  a 
certain  aggregate  inventory  I^  ,the  value  of  ^  can  be  found  from 
the  graph.  This  riiethod  is  rather  inaccurate  and  can  lead  to  very 
large  errors  in  the  estimated  ^  for  small       values  of  1°  as 
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the  curve  goes  aaymntotically  to  0, 


Fig.l 
Another  general  solution  is  by  linear  approximation 


/Vp.193  . 


*   1,0   1^  "iQi 


^=^o.AO(I*.^l2.^g^Si) 


(5a) 


where: 


-we  choose  ^°,  S?  as  an  expansion  point, 

vO      tO 


-compute  Q^,  I^  (for  given  A  .  S^)  from  equations  (3), 

(4)  respectively, 
-I^  is  the  constraining  aggregate  inventory,  and  S.  the 

forecasted  dem'and, 
-and 


A°  =  2 


u^Qf  n  -1 


Lra  4CpiS9 
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In  general,  the  aoproximation  of  equation  (5a)  could 
be  improved  by  considering  the  differential  change  in  square  root 
of  the  salee  rates  rather  than  in  sales  rates.  In  thie  case 


;^°  -^  A°  (  i:  - 


i 


^i 


(5b) 


A  better  estimate  however  is  obtained  if  solution  A  of 
(5a)  or  (5b)  is  further  refined  by  a  logarithmic  approximation; 


.2. 


A'  «  A^(l-e 
m 


An 


($) 


where  ^„  «  mixv-—- 
m       Uj 


If,  for  instance  A°  =  0,  then  Fig. 2  shows  how  the  line- 
arly  approximated  A  and  the  logarithmically  approximated  A  are 
displayed: 


1^ 


A  =  true  A  for  given  Iq 


InW   eq.    (4) 


jOg.    approx. 
eq.(6) 

.Linear  approx, 
eq.    (5a) 


Fig. 2 

It  is  obvious  that  both  solutions  by  linear  and  logarith- 
mic approximation  are  sensitive  to  the  choice  of  the  starting 
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Doint  (  A°,  S?).  The  authors  warn  about  this  /5/P.194:  "A  conveni- 
ent point  of  expansion  is  that  characterized  by  the  average  sales 
rates  and  ^°=0. . .  .  However,  this  is  not  always  an  appropriate 
point  of  expansion".  No  solution  is  offered  as  to  how  to  choose  J( 
when  A'^sQ  is  not  appropriate.  Thus,  although  the  final  formulas 
(5a),  (5b),  (6)  are  simple,  the  choice  of  A°  may  require  itera- 
tive extensive  analysis  in  order  to  obtain  a  fairly  good  estimate 
of  A*. 

It  is  apparent  that  a  more  accurate  and  efficient  method 
is  required  for  the  solution  of  (4).  The  granh  of  Iq(A)  is  shown 
in  Fig.l.  Given  a  certain  aggregate  inventory  Iq  the  problem  is 
that  of  finding  the  root  of  an  equation.  However,  it  is  easy  to 
transform  this  root-finding  problem  into  a  maximization/minimization 
which  has  the  advantage  that  some  of  the  many  optimum  seeking 
methods  available  in   the  literature  could  be  applyed  to  obtain 
the  solution  A  . 

Let  : 


I*  (7) 


(8) 


It  is  obvious  that  a  solution  to  (4)  is  obtained  when 
VJ/(A)=0,  which  is  equivalent  to  minimizing  f (  A  ) .  Fig. 3  shows  the 
relationship  among  Iq(A) , ^  (A) ,  and  f  (A) . 

Thusour  problem  is  now  one  of  unconstrained  optimization 
(minimization) , 


iM'i:=^^' 


a) 


b) 


c) 


Fig. 5 
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3.  METHODS  FUR  UNCONSTRAINED  OPTIMIZATION 

-  A  REVIEW  - 
3.1.  Several  variables  problems 

The  literature  oresents  a  large  set  of  optimum  seeking 
methods  for  multivariate  problems  /lO/,  /ll/,  /6/,  /7/. 

Fletcher  in  /6/  groups  these  methods  in  3  broad  catego- 
ries, according  to  their  basic  principle  :  gradient  methods, 
direct  senrch,and  sums  of  squares.  The  problem  to  be  considered 
is  that  of  finding  a  local  minimum  (or  if  the  function  is  uni- 
modal  the  global  minimum)  of  a  function  f(x)  where  x  is  a  vector 
of  n  variables  x=(x, jXp, . . . ,x  ) .  Let  g(x)  be  the  gradient  of  f 
with  i-th  element  9f/3x^,  and  G  the  matrix  of  second  derivatives 
with  i,  j  entry  ^^f/dx^Xy 

a)  Gradient  methods 

-  steepest  descents  in  which  the  direction  of  search 
is  s=-g.  In  practice  the  method  improves  f(x)  rapidly  on  the  first 
few  iterations  and  then  gives  rise  to  oscillatory  progress  and 
becomes  unsatisfactory. 

-  Newton'  e  method  (or  Taylor  series  method)  assumes 
that  the  function  may  be  approximated  locally  by  its  Taylor  series 
up  to  the  quadratic  terms.  Hence,  the  properties  of  quadratic 
functions  are  used  directly  to  generate  a  direction  of  search 
8=-0~  g.  The  convergence  is  rapid  if  f(x)  is  adequately  represented 
by  the  first  two  terms  of  its  Taylor  series.  However,  the  method 
has  various  disadvantages  : 

-  it  requires  the  matrix  of  second 
derivatives  to  be  provided  and  calculated  at  each  iteration, 

-  it  requires  the  solution  of  a  set 
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of  linear  equationsto  determine  each  direction  of  search  . 

Let's  note  that  when  G  is  locally  singular  the  iteration 
breaks  down;  also  if  s  turns  out  to  be  orthogonal  to  g  no  further 
progress  is  made. 

-  method  of  conjugate  directions  allows  us  to  avoid 

calculating  G;  the  cost  we  have  to  pay  for  this  is  that  a  larger 
number  of  iterations  will  be  required.  Since  G~  is  no  longer 
required  the  iteration  cannot  break  down  because  of  the  singularity 
of  0;  also,  the  method  ensures  that  the  directions  of  search  are 
downhill.  If  the  vectors  s,,82i«».»s  have  the  property 

s'^GSj  «  0   (i^4)  ;   s'^Qsj  /  0   (i=j) 
with  regard  to  a  positive  definite  matrix  G,  then  they  are  said  to  be 
conjugate.  There  are  two  ways  in  which  a  method  can  be  made  to 
generate  conjugate  directions:  the  parallel  subspace  method  and 
the  projection  method. 

All  these  methods  find  the  minimum  of  a  quadratic  func- 
tion in  a  given  number  of  searches.  For  nonquadratic  functions 
the  methods  can  be  applied  iteratively.  However,  for  nonquadratic 
functions  a  superior  method  has  been  developed  by  Davidon,  called 
the  variable  metric  method  for  minimization,  in  which  a  positive 
definite  approximation  H  to  G"  is  updated  at  each  iteration, 
and  is  used  to  generate  directions  of  search  s=-Hg. 
b)  Direct  Search  Methods  (/6/p.7,  /7/p.26) 
This  set  of  methods  is  applied  to  functions  whose  deri- 
vatives are  not  available. 

-  alternating  variable  method  -  each  variable  is 
chosen  in  turn,  all  the  others  are  kept  constant  and  the  extremum 
is  obtained  by  one  of  the  single  variable  search  methods. 
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The  method  is  very  slow,  highly  oscillatory  and  usually  fails  to 
converge. 

-  Roaenbroclc'e  modification  (or  pattern  search  as 
called  in  /8/p.B-347)  is  one  of  the  most  robust  methods  available 
for  ontimization  when  the  derivatives  are  not  available.  This 
procedure  has  been  obtained  by  imoosing  two  modifications  to  the 
alternating  variable  method  : 

1)  the  first  is  to  avoid  the  single  variable 
ootimization  for  each  direction  in  turn.  Instead  a  step  of  ore- 
determined  length  is  taken  in  each  direction  and  these  step  lengths 
are  modified  after  each  calculation; 

2)  the  second  modification  is  to  recognize 
that  the  alternating  variable  method  takes  a  large  number  of  very 
small  steps  and  then  to  try  to  avoid  this  by  realigning  the  axes. 
The  axes  are  reoriented  so  that  the  first  axis  is  along  the  most 
successful  overall  direction,  the  second  axis  along  the  next  most 
successful  direction  and  so  on.  The  change  of  axes  is  performed 
by  the  well  known  Gram-Schmidt  orthogonalization  process. 

-  simplex  method  -  the  first  step  is  to  set  ud  a  regular 
"simplex"  in  an  n-dimensional  space,  that  is  (n+1)  points  all 
equidistant  from  each  other.  The  function  is  evaluated  at  those 
Doints  and  then  the  simplex  set  is  altered  systematically  -  dropping 
some  points  and  adding  others  -  until  the  region  of  the  minimum 

is  reached.  Its  precise  location  is  found  by  interpolating  a 
quadratic  function  at  suitably  chosen  ooints.  However,  if  the 
number  of  variables  becomes  large  the  method  does  not  work  so  well. 

It  is  interesting  to  point  out  the  result  obtained  by 
Taubert  /8/  in  applying  three  search  procedures,  from  the  ones 
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enumerated  above,  for  deriving  decision  rules  for  the  aggregate 
scheduling  problem.  A  20  dimension  response  surface  was  searched 
for  the  minimum  using:  the  method  of  conjugate  direction,  Davidon's 
variable  metric  method  search,   and  the  method  of  pattern  search. 
/8/  reports  that  the  pattern  search  exhibited  the  fastest  conver- 
gence, while  conjugate  direction  search  was  the  slowest; 
Davidon's  method  yielded  an  average  performance. 

c)  Sums  of  Squares  -  here  the  special  case  is  considered 
in  which  f(x)  is  the  sum  of  squares  of  m  nonlinear  functions 
g(x) .  The  problem  can  be  solved  by  minimizing  f(x)  with  one  of  the 
methods  shown  above.  While  this  is  usually  the  safest  line,  often 
much  more  rapid  convergence  can  be  obtained  by  taking  into 
account  the  special  nature  of  f. 

When  m=n  the  problem  of  minimizing  f(x)  is  equivalent 
to  that  of  solving  a  system  of  nonlinear  equations  g(x)=0. 
This  is  an  interesting  fact  because,  given  a  system  of  nonlinear 
equations  gAjO'^O,    i«l,2, . .  ,m,to  solve  it  is  equivalent  to 
minimizing 

f(x)  *  ^   [Si^i)]^ 

If  no  exact  solution  for  minimizing  f(x)  is  available,  the 
response  surface  f(x)  may  be  searched  for  its  minimum  and  at 
least  an  approximate  solution  to  the  system  of  nonlinear  equations 
can  be  found. 

The  above  preseatntion  gives  a  general,  although  not 
exhaustive  by  any  means,  view  of  the  methods  of  searching  for  the 
minimum/maximum  of  functi  :)ns  of  several  variables.  As  a  matter 
of  fact,  there  are  almost  as  many  methods  as  there  are  researcliers 
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in  the  field. 

It  is  generally  felt  that  dimensionality  is  probably 
the  limiting  factor  in  all  these  techniques.  However,  good  nrof:ress 
has  been  made  so  far.  For  instance,  /9/ch.7  reports  a  successt\il 
search  conducted  for  finding  theminimum  of  all4-dimensional 
resDonse  surface.  The  cost  was  moderate  (about  12  minutes  on 
IBM  360/91)  and  the  authors  state:  "we  assume  that  we  have  not 
yet  reached  dimensionality  limits  so  that  the  number  of  decision 
variables  available  is  probably  somewhat  greater". 

For  the  case  of  functions  of  one  variable  many  of  the 
above  techniques  can  be  applied. However,  there  is  a  soecial  class 
of  search  methods  developed  for  single  variable  problems,  which 
are  simpler  than  the  previous  one,  and  for  particular  cases 
(ouch  as  unimodal  functions)  probably  more  efficient. 

3.2.  Single  variable  problems 

This  section  is  meant  to  give  an  indication  of  the  main 
methods  used,  with  emohasis  on  the  techniques  that  might  be 
considered  for  solving  the  HMMS  model. 

3.2.1.  Bracketing  /7/ 

It  is  most  important  as  a  first  step  in  a  optimization 
to  get  a  rough  idee  of  where  to  look  for  an  extremum;  an  useful 
idea  is  to  find  two  values  that  bracket  the  extremum.  Suppose  the 
minimum  of  a  function  f(x)  is  sought  and  it  is  known  that  the 
minimum  is  located  in  the  region  x^a.  Choose  an  increment  £  and 
evaluate  f  at  points  x,=a,  X2=x,+j^,  x,=X2'«'2^,  x.=x-+4£,  Xc=x.-»-8i, 
...,  that  is  doubling  the  increment  at  each  stage. 
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Fig.  4 

The  evaluation  is  stopned  if  either  the  minimum  is  bra- 
cketed (Fig, 4)  or  if  Xj  >X,  where  X  is  a  suitable  large  constant 
chosen  at  the  start  of  the  calculation.  The  minimum  is  bracketed 
if  at  some  stage  f(x^)<  f(Xj^_^)  and  f (x^)<  f (x^^-j^)  ;  the  bracket 
is  (Xi_n  »Xj^^]^) .  Of  course,  if  we  cannot  decide  the  search  direction 
at  the  start^both  directions  leaving  x=a  must  be  tried. 

Tf  the  first  derivative  is  available  then  if  f'(x.)<0 

J 

and  f*(x.  ^)>0  the  bracket  for  the  minimum  is  (x.,x.^,). 

If  the  value  of  f(x)  decreases  until  X  is  reached  it 
is  usually  assumed  that  the  function  is  unbounded;  if  on  the  other 
hand  the  value  of  X  is  reached  and  the  function  is  still 
increasing  then  f(a)  is  usually  taken  as  the  minimum  value. 

It  was  found  that  the  method  works  well. 


3.?. 2.  Polynomial  approximation  /?/ 

Once  a  bracket  has  been  obtained  for  the  extremum  it  is 
then  required  to  obtain  the  extremum  to  any  soecified  accuracy.  One 
simple  way  of  doing  this  is  to  use  the  information  obtained  by  the 
bracketing  procedure  directly  and  approximate  this  information  by 
a  polynomial. 
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Suppose  for  instance  that  after  a  bracketing  procedure 
developed  as  explained  in  3.2.1,  we  stopped  with  f (x^)<  f (x^^-^) , 
f(x^  , ).  For  simplicity  let  z,=x^_^,  22=x^,  Z3=x^^^.  A  quadratic 

approximation  for  f  can  be  written 

p 
f(x)  =  ax  -f  bx  ■••  c 

Using  the  known  information  the  result  can  be  summarized 

in  the  matrix  equation 

*  2 

Zn 


f(Zl) 


- 

- 

1 
1 

a 
b 

1 

c 

. 

. 

which  can  be  solved  for  a,  b,  c.  The  minimum  for  the  anoroximating 
quadratic  is  at  f'(x)=0.  Thus  an  estimate  f*af(x*)  for  the 
extremum  is  now  available.  We  can  reiterate  by  identifying  the 
best  bracket  from  the  values  z,,  Zp,  x  ,  z_,  etc. 

If  the  derivative  is  also  available  immediate  advantage 
is  obtained.  Suppose  in  3.2.1.  we  stopped  with  f'(x)<0, 
f'(Xp)>0,  hence  the  bracketing  interval  (x,  ,Xp).  Since  there  are 
now  four  nieces  of  information:  f,=f(x^),  fj^=f'(x,),  f2=f(x2), 
fl=f'(Xp)  it  is  possible  to  approximate  the  function  in  the  inter- 
val by  a  cubic: 

f(x)  =  ax^  ■♦•  bx  +  ex  +  d 
f '(x)*  3ax  +  2bx  +  c 

The  coefficients,  are  found  by  solving  for  a,  b,  c,  d 
the  following  matrix  equation  : 


■        - 

f 

a 

_ 

^1 

2 
*i    ^1 

^1 

1 

a 

'2 

"^       2 

Xp      x^ 

X2 

1 

b 

n 

3x'^   2x-^ 

1 

0 

c 

^ 

3x'J    2x2 

1 

0 

d 

^ 

« 

L 

-  17  - 

The  minimum  of  the  aDoroximating  cubic  is  x  ,  given  by 
f*(x)=0.  The  gradients  f^,  f^j,,  fi  are  compared  and  two  new  points 
are  chosen  to  give  gradients  of  opposite  sign.  The  same  procedure 
can  than  be  repeated  until  sufficient  accuracy  is  obtained. 

This  method  works  extremely  well  for  most  functions. 
However,  some  difficulty  may  arise  if  the  function  has  a  sharp 
peak. 

3.2.3.  Bolzano's  root-finding  method  /ll/  evaluates 
the  function  each  time  in  the  center  of  the  remaining  interval 
and  eliminates  half  of  the  interval  (whether  the  left  half  or  the 
right  half  is  eliminated  depends  on  the  evaluation  outcome). 
After  N  evaluations  of  this  sort,  the  ratio  of  initial  L  to  final 

0 

L«  interval  is 

^  =  2^  (9) 

Thus  the  number  of  observations  needed  to  achieve  a  given  reduction 
is  evidently 

N  =  3.321ogY^  (10) 

This  method  is  a  contender  for  solving  4'(A)=0 
defined  by  equation  (7)  and  shown  in  Fig. 3. 

Bolzano  search  /ll/requires  the  evaluation  of  both  the 
function  and  its  first  derivative.  Each  time  a  point  is  placed 
in  the  center  of  the  remaining  interval.  Let  for  instance  the 

X,+Xp 

bracketing  interval  be  (x,,Xp)  and  place  point  x-=— =^ — . 
Evaluate  f^^f  (x^^)  ,  f2=f(x2),  t^=f{x^),    f{=f(x^),    f  ^=f  (x^) ,  f^=f(x^l 
We  had  initially  f4<0,  fA>0  (since  we  are  minimizing).  The  remai- 
ning interval  is  (x^,x-)  if  f4>0,  and  (x-jXg)  if  f4<0.  Stop  when 
the  interval  is  small  enough  or  when  no  further  significant  decreasa 
in  tha  objective  functior  is  achieved. 
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3.2.4.  Direct  search  /lO/,  /7/. 

This  class  of  methods  is  concerned  with  ODtimizing  when  the 
derivative  is  not  known  or  is  comnlice ted/ inconvenient  to  be  used. 
The  idea  is  that  once  a  bracket  has  been  obtained  the  aim  is  to 
progressively  reduce  the  length  of  the  bracket  until  it  is  less 
than  a  nrescribed  limit,  or  until  no  significant  improvement  in 
the  response  function  can  be  achieved. 

Dichotomoua  search  /lO/  is  similar  to  Bolzano  search 

but  it  does  not  require  the  evaluation  of  the  first  derivative. 

Instead,  two  points  rather  then  one  are  placed  at  a  distance  £^  , 

symmetric  to  the  center  of  the  remaining  interval, 
Lo 


U/2 


tti 


«♦ 


& 


Wa   ^ 


Fig. 5 

The  evaluation  of  the  function  f(x)  at  each  point  is  called  an 
experiment.  In  Fig. 5,  after  point  x™,  x.  have  been  determined  we 

compute  f,=f(x,)  and  f>=f(x.).  Suppose  f . <  f,;  then  we  know  that 
the  minimum  of  f(x)  lies  somewhere  in  (x,,Xp),  hence  L,  is  the 
remaining  interval  after  the  first  set  of  experiments. 
After  N  experiments  (N  mu.;t  be  even  of  course)  we  can  locate  the 
minimum  within  an  interva  .  of  length: 
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h   =  '^o[^   *  ^  1  -  772  >^] 


(11) 


The  interval  S   should  be  as  small  as  possible;  it  is 
bounded  from     below  by  the  requirement  that  two  outcomes  be 
distinguishable.  It  is  important  to  point  out  that  although  the 
resolution  is  negligible  compared  to  the  original  interval  of 
uncertainity  L  ,  it  is  often  a  large  fraction  of  the  final  interval 
Lj,  if  the  search  is  at  all  efficient. 

Golden  section  and  Fibonacci  search  /7/  are  more 
powerful  than  the  dichotomous  search  technique. 

SuDoose  (a, .ap)  brackets  a  required  minimum  of  the 
function  f(x).  The  points  a,,  a>  are  symmetrically  placed  in  this 
interval,  so  that 


a,  ~  ^T ®i  ■♦■  (  1  ~^1^®2 
a^  =  (  1  -  oC^)^-^   "^  ^1®2 
and  this  division  is  illustrated  in  Fig. 6. 


i<oC<i         (12) 


Suppose  now  that  f(^^)<f(a,);  in  this  case  (a,,ap) 
brackets  the  minimum. 

Let's  take  now  .he  remaining  reduced  interval  (a,,ap) 
and  divide  it  again.  Sine  ^  the  number  of  functions  evaluations 
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must  be  reduced  to  a  minimum  it  would  be  very  convenient  to  use 
the  remaining  point  &,    in  a  further  symmetrical  division  of  the 
reduced  interval. 


<<z(^z'^i) 


<^2^*2-*3) 


Fig.  7 

Indeed,  piece  a  new  point  ae  symmetrical  to  a.  in  the  interval 
(a^.ag)  Fig. 7. 


Sc  =  (  1  -  aCo   ^^-K  ■*■   oC^a 


2=2 


"5   '  "   *^2  '°!5 
"4  ~  °^2®3  -^  (  1  -  0C2  )a 


(13) 


Since  a.  is  the  seme  in  (12)  and  in  (13)  it  follows  that 
(  1  -  cC^  )a^  +  ^i®2  ~  ^2^3  ■♦■  (  1  -  0^2   ^®? 


which  yields 


^2  = 


1«2 

1  -  eCi 


(14) 


The  method  can  be  continued  in  precisely  the  same  way 
successive  symmetric  divisions  being  performed  until  the  length 
of  the  interval  is  less  tlian  the  required  tolerance.  The  sequence 
of  fraction^  oC-,,  oCpf"    satisfy  the  recurrence  relation 


N-^: 


'N 


(15) 


The  basic  choic  is  how  to  satisfy  (15)  in  the  most  con- 
venient manner. 
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The  j^olden  section  sets  oC  =  oC,  =  orp=of„  =  , . ,  and  solving 
(15)  gives 

oC  =  0.6180335  (16) 

The  Fibonacci  search  technique  uses  the  Fibonacci 
numbers  and  works  with  a  prespecified  number  N  of  interval  divisions 
(i.e.  experiments).  Fibonacci  numbers  satisfy  : 

Fq  =  F^  =  1,   Fj^  =  F^_^   +  Fjj.2.   N^  2    (17) 

Some  of  these  numbers  are  given  in  the  table  of  Fig. 8 
for  later  use. 


N 

0 

1 

2 

? 

4 

5 

6 

7 

8 

9 

10 

11 

Fi 

1 

1 

2 

3 

5 

8 

13 

21 

34 

55 

89 

144 

Taking 


Fig.  8 


it  is  easily  seen  that  the  recursive  relation  (15)  is  satisfied. 
The  sequence  of  oCj's  for  this  method  is  : 


W     ^N-'2 

'^  ^  *N-1  (19) 

*N-1=  fj 

F 
Tn  the  final  division  with  ^v_t  =  w^  =  i  it  is  seen 

that  the  two  end  points  and  the  midpoint  have  been  selected. 

These  three  values  have  already  been  calculated  and  it  is  usual 

to  evaluate  the  function  "(x)  at  a  point  close  to  the  midpoint 

to  decide  which  half  to  clioose  as  the  final  interval. 
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A  slightly  modified  version  of  Fibonacci  search,  where 
the  last  experiment  consists  of  placing  two  points  at  a  very 
small  distance  £  from  each  other  rather  than  only  one  in  the 
middle  of  the  interval, is  developed  in  /lO/p.24. 

It  can  be  shown  that,  if  one  starts  with  a  bracketing 
interval  of  length  L  and  the  number  N  of  experiments  is  suffici- 
ently large,  the  ratio  of  the  reduced  intervalsobtained  after 
(N-1)  experiments  with  golden  section  and  Fibonacci  search  respec- 
tively is  : 


golden  section  length   ,  ,  ^^^c* 
Fibonacci  length     -i-.-l/uo 


Thus  for  large  N  the  Fibonacci  search  gives  a  17%  better 

result  than  the  golden  section.  Both  methods  are  superior  to  the 

dichotomous  search  because  at  every  step  the  ooint  remaining  in 

the  reduced  interval  is  used  in  a  further  symmetrical  division; 

this  feature  If^ads  to  a  better  use  of  the  information  available 

after  each  experiment. 

Indeed,  after  N  experiments  the  initial  interval  L  is 

o 

reduced  down  to: 


^N  '^l*^2---^N-l^o 
For  the  golden  section 


Lh  = 


(1.618)^"^ 


and  with  Fibonacci  search 

If  these  results  are  compared  with  the  result  of  the 
dichotomous  search  (equation  11)  it's  obvious  that  the  golden 
section  and  Fibonacci  search  oerform  much  more  efficiently. 

For  an  unimodal  function  it  can  be  oroven  that  both 
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golden  section  and  Fibonacci  search  technique  always  work. 

Lattice  search  by  Fibonacci  technique  /lO/ 

Wilde  /lO/  raises  the  point  that  it  may  sometimes  be 
advantageous  to  convert  an  ordinary  continuous  search  into  a  lattice 
search  artificially. Indeed,  the  result  of  a  continuous  Fibonacci 
search  willbe  an  interval  that  contains  the  optimum.  When  one  is 
expected  to  make  a  decision  based  on  the  results  of  a  search  it 
is  a  bit  frustrating  to  be  confronted  with  an  interval  of  uncer- 
tain ty.  A  oreciee  ooint  would  be  oreferable,  since  a  specific 
decision  is  called  for.  One  could,  of  course,  choose  a  ooint  at 
random  in  the  final  interval  of  uncertainity,  but  most  people 
would  prefer  a  point  where  a  measurement  had  already  been  made. 
Thus,  to  avoid  these  difficulties,  the  original  problem  can  be 
converted  into  a  lattice  problem  by  placing  a  number  of  points  in 
the  bracketing  interval  so  that  the  final  answer  will  be  a 
enecific  point  on  which  a  firm  decision  can  be  based. 

Suppose  it  has  been  decided  that  a  number  N  of  exneri- 
mente  will  be  performed  to  search  by  Fibonacci  a  bracketing 
interval.  Partition  the  interval  into  F«  units  using  F«-l  points, 
not  necessarilly  equidistant.  These  points  form  a  lattice.  Let's 
associate  the  lattice  points  with  the  integers  1  through  Fjr-1. 
Thus,  we  are  dealing  here  with  an  original  interval  of  length  F^, 

units.  According  to  relation  (19)  and  Fig. 6  the  first  two  exoeri- 
ments  should  be  placed  at  a  distance  of  ^-iFm^Fm-,  so  that  the 
first  two  experiments  will  coincide  with  points  of  the  lattice, 
namely  ooint  number  F«  ^  and  point  number  Fjt-F«t_-j=Fm  2* 

Since  the  length  of  the  interval  remaining  is  also  a 
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Fibonacci  number,  we  see  that  the  third  experiment  will  also  fall 
on  one  of  the  lattice  points.  This  procedure  may  be  continued 
until  N-1  experiments  have  been  used  up  and  the  length  of  the 
interval  of     uncertainity  is  dawn  to  oC^cCg. . .  fl^jij_2^N~^2"^  units. 
The  sequence  is  stopped  and  the  unique  point  left  inside  the  two 
unit  long  interval  is  compared  to  the  end  points  and  the  best  will  be 
the  estimate  of  the  optimum,  on  which  then  all  decisions  will  be  based 
Let  us  emphasise  again  that,  in  order  for  this  tec^hni- 

que  \Q  woric  properlvi  the  JLnilial  bracltgting  interval  must  be 

partitioned  bv  a  number  of  points  equal  to  a  Fibonacci  number  leas 

4.  SOLUTION  TO  HMMS  MODEL  BY  LATTICi;;  SEARCH 
WITH  FIBONACCI  TECHNIQUE 

From  all  the  search  methods  reviewed  in  section  "5 
Kibonacci  technique  gurantees  the  largest  interval  reduction  in  A 
given  number  of  steps.  Moreover  it  does  not  require  the  use  of 
derivatives,  a  fact  that  is  advantageous  in  our  case  because  the 
repeated  evaluation  of  the  derivative  of  (8)  would  require 
additional  computing  effort  given  the  fact  that  it  is  computatio- 
nally more  complicated  than  the  original  function. 

Thus,  the  task  of  this  section  is  the  search  for  the 
minimum  of  function  (8)  whose  graph  and  equation  will  be  repeated 
here  for  convenience  (Fig. 8). 

We  will  use  the  same  example  as  the  one  given  by  the 
authors  in  /Vp.l96  in  order  to  be  able  to  compare  performance. 
Assume  that  the  estimates  of  the  setup  costs  and  the  costs  of 
holding  inventory  for  each  product  are  reviewed  and  revised 
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annually.  Forecasts  are  also  made  of  the  average  monthly  sales 
rate,  S' ,  for  each  oroduct  for  the  coming  year.  The  aggregate 
inventory  levels  are  planned  with  view  to  both  labor 
requirements  and  costs  associated  with  inventory.  The  aggregate 
inventory  will  therefore  be  measured  in  labor  hours  by  multiolying 
the  units  of  each  product  by  the  conversion  factor  u^,  which  has 
the  dimension  labor  hours  per  item.  The  relevant  data  are  summa- 
rized in  the  table  of  Fig.  9. 


minimize  f(^) 


f(A)  = 


>».  ^ 


Fig.  8 


Product 

Units 

1 

2 

3 

"li 

1 

1 

2 

;S/ it  em-month 

°pi 

10 

10 

30 

SJ/lot 

"i 

1 

5 

2 

hours/ item 

Si 

2000 

2000 

10000 

items/month 

(8) 


Pig. 9 


The  point  of  expansion  for  linear/logarithmic 
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\0. 


aporoximation  was  chosen  ^  =0  and  the  averages  sales  rates  for 
each  item.  The  table  in  Pig. 10  shows  alternative  estimates  of  ^ 
and  the  error  in  Iq  compared  to  the  imposed  aggregate  level  Iq 
(Iq  is  computed  using  the  estimate  of  A  ). 


r;. 

800 

1500 

Is 
nttutHnq 

£tror 

o/ 
/o 

^  by  linear 
aDproximation  (5a) 

-0.22 

931 

16.4 

0.22 

« 

« 

A  by  logarithmic 
aporoximation  (6) 

-0.41 

832 

4 

0.14 

1611 

7.4 

>i*by  graphical 
method 

-0.45 

815 

1.9 

0.12 

1481 

-1.3 

infeasible  since  the  estimate 

•Y   .  ^li 
of  .22  exceeds  ^_=min-i-=  .20 
m     u> 


Pig. 10 

Both  linear  and  logarithmic  approximations  overstate 
>*  and  may  lead  to  considerable  error,  even  to  infeasibility. 
The  graphical  estimate  is  more  accurate  but  it  does  not  seem  to 
be  very  appealing  from  a  practical  ooint  of  view  because  in  order 
to  get  a  fairly  good  curve  one  needs  a  large  number  of  ooints 
to  be  evaluated;  moreover,  as  shown  in  section  2,  given  the  asymp- 
totic characteristics  of  f ( ^ )  the  errors  in  reading   off  the 

m 
graph  for  small  values  of  I-  are  likely  to  be  quite  large. 

Let  us  now  estimate  ^  using  search  techniques. 

4.1.  Bracketing  A* 
4.1.1.  The  case  I*<Iq(0) 

In  this  case  ^  <  0.  An  upper  bound  on  ^  is  easily 
obtained  if  we  observe  th ^t  the  intersection  between  the  tangent 
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to  y^CA  )  at  ^=0  and  the  horizontal  axis  is  always  larger  than 
y* ,   Let's  call  this  upper  bound  ^  .  The  reason  behind  the  fact 
that  we  used  ^(^  )  rather  than  the  searched  response  function 
f ( A  )  to  derive  bounds  will  become  clear  in  section  5.  Anyhow, 
the  brackets  developed  with  *^  (  ^ )  are  perfectly  valid  for  f(^  ), 
given  the  equivalence  shown  in  Fig. 3. 


^li-  ^i  "  Q 


(7) 


Fig. 11 

A  lower  bound  can  be  obtained  using  the  bracketing 
technique  described  in  3.2,1.  Use  i~-^     and  go  out  f rom -^ 
with  a  double  increment  2^  (see  Fig. 11).  Evaluate  U'C  A^-2^)  and 
repeat  the  procedure  (if  necessary)  with  4^,  8^  etc.  until^^i^) 
becomes  negative.  Then  stop  with  the  lower  bound  Ap, 

4.1.2.  The  case  I*>Iq(0) 

In  this  case  0  <  ^*<  X^,   The  interval  (0,  Ti^)    is  a 
possible  choice  for  bracketing  ^  .  However,  for  the  porpose  of 
computing  an  upper  bound  on  the  error  of  estimation  (see  section  6), 
this  interval  is  unsuitable.  An  uoper  bound  better  than  ^^   is 
required.  Two  subcases  can  be  distinguished  here  : 

a)  A  corresponi  ing  to  the  intersection  between  the 
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tangent  to  j(.A)    at  J\=0  and  the  horizontal  axis  is  smaller  than 
/(_  (Fig. 12).  In  this  situation  the  intersection  is  labeled  A  , 
and  a  convenient  bracketing  interval  is  (0,  ^,,) 


Fig. 12 

b)  A  corresponding  to  the  intersection  between  the 
tangent  tp  U/(^)  at  A  =0  and  the  horizontal  axis  is  larger  than 


>^  (Fig. 15). 


Fig. 13 

Then,  a  good  bracketing  procedure  is  a  binary  search  /12/p.82, 
The  algorithm  is  the  following  : 


Let  k=l,  X=o 

'     0 

Step  1  -  Evaluate  U^i  :X  )    et 


^m"  ^k-1 


If  (Pi  ^^)  ^O  go  to  step  II. 
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If  ^(  ^^)=0   BtoD  with  5^*=  >^. 

If  a/(  \)>0   stop  with  ;^^=  ^  and  \=  \_i 

Step  2  -  Let  lc=k+l.  Go  to  step  I, 

In  Fig,  13  the  algorithm  terminated  after  two  evaluations, 
yielding  (  A,,  ^  )  as  the  bracketing  interval. 

Let's  note  that  the  binary  search  is  actually  the 
bracketing  technique  of  section  3.2.1.  applyed  in  reversed  order; 
this  was  Doesible  because  we  started  out  with  a  finite  interval 

(0,  AJ. 

4.2.  Estimating  ^ 

Once  the  bracketing  interval  has  been  obtained, 
Fibonacci  search  technique  will  be  applyed  to  an  artifficially 
constructed  lattice.  Namely  we  decide  first  to  perform  N 
Fibonacci  exoeriments;  consequently  F„,  points  will  partition 
the  initial  interval  into  F„  units,  and  the  technique  is  apolyed 
to  this  lattice.  The  issue  of  how  to  choose  N  will  be  addressed 
in  section  6. 

4.3.  Numerical  examples 

Example  1 

Consider  the  example  solved  by  the  authors  for  an 
aggregate  inventory  of  Iq=800.  All  relevant  data  are  given  in  the 
table  of  Fig.  9 

In  our  perticulnr  case  we  have: 

iU(  1  >i    -  1  iMOOOO"  ^  5  lAouOO  ^  1/600000   ann 
The  tangent  to  ^(  J^  )  at  A  =0  is: 
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0/-  348  =  1573.86  ^ 
Applying  the  bracketing  procedure  shown  in  4.1.1.  one 
obtaines: 

:X^   =  -.22   ;  }i^   «  -.66 

hence  the  inteval  within  which  to  search  for  ^^  is  (-.66, -.22)^ 
the  response  function  whose  minimum  we  are  searching  for,  is  : 

f(^)  ^\f{^   )| 
and  it  was  pictured  in  Fig,  8. 

Choose,  for  instance,  to  perform  N=5  Fibonacci  experi- 
ments. Consequently,  the  interval  (-.66, -.22)  will  be  partitioned 
into  Fjj=8  units  using  Fj^_-,=7  points  equally  spaced.  However,  there 
is  no  special  requirement  that  the  points  be  equidistant;  con- 
sequently, if  we  have  any  suspicion  that  A  might  lie  in  a 
certain  subinterval  of  the  bracketing  interval  we  might  want  to 
distribute  the  nointe  closer  in  that  subinterval  and  further 
apart  in  the  rest. 

The  increment  of  the  lattice  is  h  =  -^4^  =.055. 


"W 


h^.OSS 


\ 


7       I      I      ff     / 

— i ^     ^ ^^ ®— 


I ^ 9 ^ ^ -^ • ^ 


-.6C    -.«W   -yS"   -.♦9y  -H   -'iW  -.5  J   -.275*  -.2^ 

Fig. 14 

First  and  second  Fibonacci  experiments 
Place  two  point;  at  a  distance  oC^Fj*  =  F«  -,  units 
from  the  ends  of  the  interval;  Fjj_,  =  5,  so  the  two  points  are 
point  no. 5  (^=-.385)  and  point  no. 3  (^=-.495). 

f(-.495)  =  2.0324,    f(-.385)  =  42.7343 
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Third  experiment 

The  remaining  interval  i8(-.66,-.385) ;  olace  in  the 
interval  a  point  symmetric  to  (-.495).  This  point  will  be  A =-.55. 
f(-.55)  =  21.53805 
Fourth  experiment 

The  remaining  interval  is  (-.55, -.385)  and  the  point 
chosen  iAa-,44 

f(-.44)  =»  19.27695 
Fifth  experiment 

The  remaining  interval  is  (-.55, -.44)  and  it  contains 
only  one  point  ^=-.495.  This  is  the  best  point  in  the  lattice,  so 
it  is  chosen  as  our  estimate  of  A  . 
A*  =  -.495 
Fig.  14  shows  the  five  Fibonacci  experiments  indicating 
the  sequence  in  which  they  have  been  chosen. 
Example  2 

Consider  the  same  data  except  for  the  aggregate  inventory 
level  which  shou]J  be  now  I*=1500. 

il/r  :A  ^  -  1  \/40000  .  5  yMOOOO  ^  \/ 600000   ,  ^^. 

The  tangent  to  ^(^  )  at  ^=0  is: 
0^  +  352  =  1573.86  A 

The  vertical  asymotote  for  ^( A )  is  ^  =0.2. 
The  intersection  of  the  above  tangent  with  the  horizontal  axis 
^-0   falls  beyond  ^^^=0.2  at  ^  =  .22,  hence  this  point  is  infeasible 
and  can't  be  used  as  an  upper  bracket  (see  Fig. 13).  Apply  then  the 
binary  search  described  in  4.1.2.b, 


k  =  1  , 


32 


>o  =  ° 


^(  \)    =  -110.1339  <  0 

V^i 


k  =  2 


^2  =  -^1  * 


=  .15 


^(  ^2)  =  202.5537  >  0 

Hence,  the  bracketing  interval  is  ( -^^ ,  ^o^'^^     &*  ^u^~ 
=(.10, .15). 

Let's  choose  to  perform  N=5  experiments.  The  corresoonding 
Partitioning  of  the  interval  is  shown  in  Fig. 15.  After  the  search 
the  best  estimate  turns  out  to  be  A  =.125.  The  Roman  numerals 
in  Pig, 15  show  the  sequence  in  which  the  experiments  have  been 
placed. 


h-.0062^ 


4 


^   I  ^  i 


V 


-#- 


4^ 


,iO        .i96i^    'Wr        >Hf75-     J2T      Jms   Ji7f       J^iff     ./5" 

Fig.    15 


In  example  1  the  estimated  ^  =-.495  yields  an  aggregate 
inventory  1^=798,  i.e.  an  error  of  0.25%  as  comoared  with  the 
imposed  Iq=800.  In  example  2  our  estimate  ^  =.125  leads  to 
Iq=1509,  i.e.  an  error  of  0.6%  relative  to  the  constraining 


i;=i5oo. 


The  point  which  we  would  like  to  emphasize  here  is  not 
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primarily  the  accuracy  itself  obtained  above;  this  accuracy  is 
more  or  less  important  if  we  look  at  it  in  connection  with  the 
accuracy  of  the  input  data  (setup  and  holding  costs,  demand 
forecasts,  conversion  factors).  The  imoortantooint  is  rather  the 
small  computational  effort  required  to  obtain  this  accuracy,  and 
the  fact  that  it  can  be  easily  improved  by  increasing  the  number 
of  exneriments  (points)  by  just  a  few.  Indeed,  an  initial  interval 
of  uncertainity  can  be  reduced  to  less  than  one  per  cent  of  its 
original  length  after  only  eleven  sequential  experiments. 

"5.  AN  ALGORITHM  FOR  FINDING  THE  ROOT  UF  AN  £QUATIUM 
BY  LATTICE  SEARCH  WITH  FIBONACCI  NUMBERS 

Fibonacci  search  technique  is  meant  for  finding  the 
minimum/maximum  of  a  function  of  a  single  variable  (unimodality 
is  desirable  to  ensure  a  successful  search)  This  is  why,  although 
our  original  problem  was  to  find  the  root  of  a  complicated  equa- 
tion (7),  we  transformed  it  to  obtain  a  peaked  function  f(^) 
(see  Fig. 3)  on  which  the  search  has  been  performed. 

In  this  section  the  original  problem  will  be  addressed. 
Wo  propose  a  general  purpose  lattice  -  Fibonacci  type  algorithm 
for  searching  for  the  root  of  an  equation  within  some  previously 
determined  bracketing  interval.  This  method  will  be  used  to  search 
the  HMMS  model  and  it  will  be  clear  that  it  requires  a  smaller 
number  of  function  evaluations;  the  number  of  evaluations  deoends 
on  where  in  the  bracketin/^  interval  is  the  root  located. 

In  Fig,  16  ,  (A,B)  is  a  bracketing  interval  for  A  ; 
y     is  at  the  same  time  the  root  of  U^(^)=0  and  the  minimum  of 
f(^).  Suppose  C,D  are  th*  first  two  lattice  -  Fibonacci  experi- 
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mente  placed  in  (A,B);  then,  according  to  the  technique  in  section 
4,  the  remaining  interval  will  be  (C,B)  as  f(D)^f(C). 


Fig. 16 

However,  if  we  look  at  ^(  ^)    instead  of  f ( ^ )  the 
remaining  interval  would  be  (C,D)  ae  Q'(C) 'I'(D)  <  0.  Hence  the  idea 
of  modifying  the  criterion  of  selecting  the  remaining  interval. 
while  keeping  the  rule  of  placing  experiments  according  to  the 
lattice  -  Fibonacci  search  method.  We  realize  that  the  root  A 
could  fall  not  only  in  the  interval  (C,D),  but  also  in  (A,C)  or 
(D,B),  as  shown  in  Fig,  17,  so  after  an  experimental  step  the 
remaining  interval  could  be  either    (A,C),  (C,D)  or  (D,B). 
In  order  to  further  apply  lattice  -  Fibonacci  technique  on 
the  remaining  interval  it  must  contain  a  number  of  lattice  points 
equal  to  a  Fibonacci  numb(>r  less  1.  Let's  investigate  (A,C), 
(C,D),  (D,B)  from  this  point  of  view. 
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Fig. 17 


PROPOSITION  Let  the  open  interval  (A,B)  contain  F^^-l  points, 
where  F„  is  a  Fibonacci  number.  Suppose  the  first 
two  experiments  C<D  have  been  placed  in  (A,B) 
according  to  the  lattice  -  Fibonacci  search 
technique.  Then,  each  open  interval  (A,C),  (D,B) 
will  contain  F  -1  points,  and  the  open  interval  (C,D), 
F^-1  points,  where  F„,  F^  are  Fibonacci  numbers. 


Proof  Points  A,  B,  C,  D,  are  displayed  like  in 

Fig  16.  Let  N^„  denote  the  number  of  points  contained 

in  some  open  interval  (iftS), 

According  to  3.2.4.  ooint  D  is  the  Fjj_,-th  Doint 

from  A,  and  C  is  the  F«  -.-th  point  from  B.  Thas 


\d  =  ^N-1  -  1 


^CB  =  %-l  '  ^  • 


Consequently  : 
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NaC  =  \b   -  ^cb  -  1  =  (F^  -  1)  -  (Fjj.;^  -  1)  -  1  = 

Similarly  : 

NCD  =  \d   -  ^AC  -  1  =  (^N-l  -  ^^  -  ^  V2  -  1)  -  1  = 

and  clearly  Nj^q  «  N.^  . 


COROLLARY  -  Let  ^i^  )    be  a  continuous  function  of  ^    with  a  root 
contained  in  some  bracketing  open  interval  (A, 3). 
Assume  (A,B)  contains  F^-1  ooints,  where  F  is  a 
Fibonacci  number,  and  that  the  first  two  experiments 
C<D  have  been  placed  in  (A,B)  according  to  the  lattice 
-  Fibonacci  search  technique. 

If  the  next  (reduced)  bracketing  interval  is  chosen 
such  that  ^f{ir)^iS)^0,     y,  cr€{A,B,C.D] 
then  Nj-^  is  equal  to  a  Fibonacci  number  less  1. 

Proof  -  The  proof  follows  directly  from  the  above 
proposition  coupled  with  figures  16,  17. 


Of  course,  from  all  possible  intervals  (0  |0  ) 
constructed  such  that  ^^  (  J^)  ^(cT  )  ^0, 

}^fO£,  ?A,B,C,DS  we  will  choose  the  one  with  the  minimum  Ny«jr  for 
the  purpose  of  efficiency. 

Summary  :  In  solving  the  HMMS  model  we  will  search  for  the  unique 
root  of  ^(^)\    a  lattice  will  be  constructed  inside  the  initial 
bracketing  interval  using  a  number  of  points  equal  to  a  Fibonacci 
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number  less  one.  At  each  step  two  experiments  are  olaced  accor- 
ding the  method  described  in  3.2.4.,  and  the  remaining  interval 
will  be  chosen  as  specified  by  the  corollary  above.  The  procedure 
will  be  then  repeated  with  the  remaining  bracketing  interval. 

Example 

For  the  purpose  of  comparison  the  examnle  already  solved 
in  4.3.  will  be  tackled  by  the  modified  algorithm.  Consider  the 
data  in  the  table  of  Fig. 9,  and  an  aggregate  inventory  of  I*=1500. 

The  bracketing  interval  is  ( .10, .15) ( see  4.3);  we  have 
already  chosen  N=5,  so  there  will  be  Fjj-1=7  points  in  (.10,. 15). 

I   jr  I 

A  CD  fi 

-U — ^. — '.     \    %    %    ''    '''     U 

.to    'WA^S"  •'««"  **^''^     '^^      ''^'^^   ''3^^  ./Vi^r  ,/y       ^ 

Fig. 18 

First  and  second  Fibonacci  exoerimenta  are  placed  at 
Fjj_i  units  from  both  ends.  So,  C*. 11875,  D=. 13125 
'I^C. 11875)  =  -25.5509 
^(.13125)  =  47.7335 
Evidently,  since  ^^fC)  0^(0)^0  the  remaining  interval 
is  (C,D). 

Third  Fibonacci  experiment  -  interval  (C,D)  only  contains 
one  point,  so 

^^(.125)  =  8.94<|y(C)|,a'(D)  . 
Hence   ^  =.125,  and  with  the  modified  method  3  function 
evaluations  were  needed  iistead  of  4  in  section  4.3.  The  modified 
technique  is  faster  in  reducing  the  interval  of  uncertainity 
especially  i^t  the  very  beginning.  This  is  easily  seen  if  we  compare 
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the  remaining  intervals  obtained  with  the  two  methods  applied 
to  the  lattice  shown  in  Fig.  18: 


Lattice-Fibonacci 


Exoeriment 
performed 

^  =  .13125 
^=.11875 


^  =  .125 


Remaining 
interval 


(.10,. 13125) 


(.11875,. 13125) 


Solution 


Modified  technique 

Experiment   Remaining 
performed    interval 


^  =  .13125 
>  =  . 11875 


/\  =  .1125     (.1125,. 13125)   /  -^  =  .125 


(.11875,. 13125) 


Solution 


The  speed  of  convergence  of  the  oropoeed  technique 
is  influenced  by  the  position  or  the  root  in  the  bracketing 
interval:  the  fastest  reduction  in  the  uncertainity  interval  is 
obtained  when  the  root  lies  in  the  central  oart  of  the  bracketing 
interval.  The  results  of  two  examoles  worked  out  on  a  lattice 
divided  into  89  intervals  (units)  are  shown  below.  In  the  first 
example  it  was  assumed  that  the  root  is  located  in  the  central 
Dart,  while  in  the  second  the  root  was  located  at  the  left  extreme. 

Example  with  root  centrally  located 


Lattice 

-Fibonacci 

Modified 

technique 

Experiment 

Remaining  inter- 

Experiment 

Remaining  inter- 

performed 

val,  in  units 

performed 

val,  in  units 

I,  II 
III 

55 
34 

I,  II 

21 

IV 
V 

21 
13 

III,  IV 

5 

VI 

VII 

VIII 

8 

5 
3 

V,  VI 
VII 

2 

Solution 

IX 

2 

X 

Solution 

-  -^9  ^ 
Example  with  root  at  the  extreme  left  end 


Lattice-Fibonacci 


Exoeriment 

Remaining  inter- 

performed 

val, 

in  units 

I,  II 

55 

III 

54 

IV 

21 

V 

13 

VI 

8 

VII 

5 

VIII 

3 

IX 

2 

X 

Solution 

Modified  technique 
Experiment   Remaining  inter- 
performed    val,  in  units 


I,  II 

III,  IV 

V,  VI 

VII,  VIII 
IX 


34 


13 


Solution 


The  advantage  of  the  fast  convergence  at  the  beginning 
of  the  search  could  be  exploited  by  checking  the  uncertainity 
interval  after  each  iteration.  We  do  not  have  to  perform  all  the 
experiments  required  by  a  Fibonacci  search;  we  can  stop  whenever 
the  uncertainity  interval  is  accurate  enough. 

Bxample  -  Consider  again  Fig. 18,  where  the  initial  bracketing 
interval  is  (.10,. 15).  We  place  2  points  in  the  interval: 
■^,=.11875,  ^2^*1^125  and  the  remaining  interval  in  the  modified 
technique  is  (.119,. 13).  Suppose  this  is  sufficiently  accurate, 
so  stop  with  ?r   (.119,. 13). 

In  order  to  obt^^in  the  same  interval  of  uncertainity 
the  original  Fibonacci  te-.hnique  would  have  required  4  function 
evaluations  at  :  ^^=.119,  ^2~*^'»  '^3=»^38,  X  =  .126. 

It  is  also  apparent  from  the  results  presented  above 
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that  a  smaller  number  of  reduced  intervals  must  be  kept  track  of, 
which  reduces  computer  time. 

Thus,  we  will  resort  to  the  modified  lattice-Fibonacci 
technique  in  the  computer  program  for  solving  the  HMMS  model, 
using  as  a  response  function  'rC^  )=0. 

6.  CHQQSINQ  THE  NUMBER  OF  POINTS  IN  THE 
BRACKETING  INTERVAL 

There  are  two  requirements  to  be  met,  and  which  place 
bounds  on  the  number  of  points  dividing  the  initial  bracketing 
interval: 

a)-  the  requirement  that  two  adjacent  function 

evaluations  be  distinguishable;  this  places  an 
upper  bound  on  the  number  of  points  in  the  interval, 
b)-  the  accuracy  of  the  result  which  is  better  if  the 
number  of  points  is  larger;  this  requirement  sets 
the  lower  bound  on  the  number  of  points  . 
Let  (A,B)  be  the  initial  bracketing  interval  and  N.^ 
the  number  of  Doints  partitioning  the  open  interval  (A,B). 
Then  /lO/p.37: 


where 


"ab  «  ^n  -  1 


N  =  the  integer  part  of  4.785  log  ^  -  0.328 
£,-   the  minimum  spacing  for  which  two 

outcomes  are  distinguishable 
F^  =  Fibonacci  number 


However,  given  the  available  modern  computing  capability  it  is 
probably  unlikely  for  thi.-,  upoer  bound  to  be  constraining  in  the 


-  41  - 
case  of  a  production  problem  like  the  HMMS  model.  In  such  cases 
we  expect  the  cost  of  the  computational  effort  to  place  an  upper 
bound  on  the  number  of  points.  Thus,  we  willtend  to  select  the 
smallest  number  of  points  allowed  by  the  accuracy  requirement. 

One  way  to  approach  this  problem  is  the  following:  given 
the  Inttice  increment  h  and  the  fact  that  our  estimate  of  }i 
is  one  of  the  lattice  points,  we  know  that  we  can  not  be  off 
the  true  root  by  more  than  h.  The  change  in  the  aggregate  inventory 
corresponding  to  a  variation  of  h  in  A  can  be  approximated  by 

^-le-  h,  for  h  small.  As  we  are  interested  in  the  proportional 
change  in  aggregate  inventory,  we  will  use  : 

1  ^^Q 
Ej-f-^h  (22) 

where  Ey  is  the  proportional  error  in  aggregate  inventory. 

Clearly,  setting  an  unper  bound  on  the  error  E^  is 
equivalent  to  limitting  h  from  above  and  N.^  from  below.  But 
Ey  is  a  function  of  ^    and  obviously  Ej(  ^  )  can  not  be  computed 
before  A  is  found.  How  should E^  be  computed  in  order  to  deter- 
mine the  minimum  number  of  points  needed  ? 


«1   V2CpiSi 

372 


diQ  _  t^  *  (Cji  -;\ui) 
^Mc,i-.:^u,)^^^ 


1     -Lq  d;\         ,ji^       V^CpiSi 


Analyzing  expression  (23)  we  find  that  Ej( A  )  increases 

Ti 
monotonically  from  0  to  e>o  ,  for  -oo<.y\^   min  .  Consequently, 

an  •      ^^ 

givenAinterval  (^pi^y)  bracketing  ^  ,  where  ^ -  ,  ^   are 
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the  lower,  and  upper  ends  of  the  interval,  respectively, 

Ej(^^)>Ej(  ;\*)   .  (24) 

Hence,  the  error  computed  at  ^   is  an  upper  bound  on 

the  error  in  the  aggregate  inventory  corresponding  to  the 

estimated  ^  . 

L 
If  L-  is  the  initial  bracketing  interval  then  h=~ 

80  taking  (22)  and  (24)  into  account  we  can  find  the  lower 

bound  on  F^jt 

where  all  elements  of  the  right  hand  side  are  computed  at  .^  = ^  , 
and  Ej(  A  )  is  specified  by  the  user  of  the  model. 

Example 

Consider  the  data  shown  in  the  table  of  Fig. 9  and  an 
aggregate  inventory  level  of  I-.  =  1500.  In  section  4.3  the 
bracketing  interval  for  ^  was  found  to  be  (.10,. 15).  In  how 
many  units  should  we  partition  this  interval  ? 

Assume,  for  instance,  that  we  want  to  limit  the  error 
in  aggregate  inventory  to  5%,  i.e.  0.05.  Consequently  we  set 
Ej(  >^)  =  Ej(.15)  =  0.05   . 

Iq(.15)  =  1702.5 

^(   =.15)  =  10556.3 

\   =  -05 
Applying  (25)  we  obtain 
Fjj^  6.083 

Since  Fj^  must  be  a  Fibonacci  number  we  look  it  uo 
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in  the  table  of  Fig. 8  and  find  Fj,  =  8  corresponding  to  5  Fibonacci 
exoeriments.  The  number  of  points  in  the  bracketing  interval  will 
be  Fj^-1  =  7. 

7.  SOLUTION  BY  BOLZANO'S  ROOT  FINDING  SEARCH  METHOD 

As  shown  in  3.2.3.  i  after  N  evaluations  of  the  response 

function,  the  initial  bracketing  interval  L  can  be  reduced  to 

LV2  with  Bolzano's  method.  If  we  have  dealt  with  a  continuous 
o 

*%  N   .     .       .  ... 

^,  after  comparing  2  with  Fibonacci  numbers  in  Fig. 8  it  would 

have  been  quite  obvious  that  Bolzano's  method  achieves  a  greater 

interval  reduction  in  a  given  number  of  iterations  N.  However, 

for  reasons  shown  on  page  23,  we  preferred  to  turn  our  problem 

into  a  discrete  version,  in  which  case  it  is  not  immediately  clear 

which  method  perfonns  bettpr. 

The  general  approach  with  Bolzano's  technique  is  similar 

to  what  we  did  so  far: 

-  set  the  maximum  acceptable  relative  error  in 
aggregate  inventory^ 

-  determine  the  number  of  intervals  in  which  the  initial 
bracketing  interval  must  be  oartitioned;  the  number  of  intervals 
should  be  equal  to  a  power  of  two,  so  the  oroblem  is  that  of 
finding  the  smallest  k  in: 

-  apply  Bolzano  search  on  the  lattice  just  constructed. 
A  computer  program  has  been  set  up  to  experiment  with 

the  two         methods:  the  modified  algorithm,  and  Bolzano's 
search.  The  following  computer  report  shows  the  result  of  an 
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example  run  with  the  data  in  the  table  of  Fig. 9,  for  an  imoosed 
aggregate  inventory  of  800. 


WHAT    IS    THE    NUMBER    OF     ITEMS? 

13 

WHAT     IS    THE    AGGREGATE    INVENTORY? 

WHAT     IS    HIGHEST    ACCEPTABLE    ERRO".    -    PERCENTAGE    ? 
13 


MODIFIED    LATTICE-FIBONACCI     SEARC?J    ALGORITHM 

NUMBER    OF    ITEMS    =       3 

IMPOSED    AGGREGATE    INVENTORY    =       800 

MAXIMUM    ADMISSIBLE    ERROR    =       3    % 

LENGTH    OF    INITIAL    BRACKETING    INTERVAL    =        .4A1869 

ACTUAL    NUMBER    OF    ITERATION?    =       4 

ITERATIONS    REQUIRED    BY    PURE    FIBONACCI     SEARCH    =       6 

ESTIMATED    LAGRANGE    MULTIPLIER    =       -.ii9P.85A 

RESULTING    AGGREGATE    INVENTORY    =       798.762 

ACTUAL    ERROR     IN    AGGREGATE    INVENTORY    =       -.15477    % 

PRODUCTION  PLAN 

ITEM       1  163  UNITS 

ITEM       2  107  UNITS 

ITEM       3  448  UtJITS 


LATTICE    SEARCH    BY    30LZArJ0    METHOD 


NUMBER    OF    ITERATION?    =       4 

ESTIMATED    LAGRANGE    :'ULTIPLIER    =       -.497103 

RESULTING    AGGREGATE     INVENTORY    =       797.189 

ACTUAL    ERROR     IN    AGC  IGATE    INVENTORY    =       -.35  13  64    % 

PRODUCTION  PLAN 

ITEM       1  1  63    U[   ITS 

ITEM       2  107    UflTS 

ITEM       3  .      447    U^ ITS 
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It's  clear  that  both  the  modified  technique  and  Bolzano 
method  had  the  same  oerformance,  while  outDerforming  the  oure 
Fibonacci  search. 

The  explanation  for  the  equal  performance  is  the  following; 
the  lower  limit  on  the  number  of  partitioning  intervals  was  a  num- 
ber less  than  13.  Consequently,  by  formula  (25)  a  number  of  13  di- 
visions has  been  considered  for  the  modified  algorithm,  while  by 
formula  (26)  a  number  2  =16  divisions  were  constructed  for  Bolzano's 
search.  Thus,  although  Bolzano  technique  is  in  general  the  fastest, 
the  discrete  nature  of  our  oroblem  lead  to  some  loss  of  efficien- 
cy. However,  for  larger  intervals,  or  for  a  lower  desired  error, 
Bolzano's  method  performs  better  than  the  modified  lattice-Fibo- 
nacci algorithm  (see  attached  computer  reports). 

Conclusion  -  For  inventory  levels  Iq  smaller  but  relatively  close 
to  Iq  (A=0),  and  for  I*   larger  than  Iq  (A=0)  the  initial  bra- 
cketing interval  is  expected  to  be  small,  and  there  is  no  signi- 
ficant performance  difference  between  the  modified  lattice-Fibo- 
nacci search  and  Bolzano's  technique.  Both  perform  extremely 
efficiently.  For  I^  considerably  smaller  than  I^  (A=o)  Bolzano's 
method  oerforms  better  and  is  to  be  preferred. 

8.  A  COMPUTER  PROGRAM  FOR  THE  MODIFIED  LATTICE-FIBONACCI 
SEARCH  TECHNIQUE  AND  FOR  BOLZANO'S  METHOD 

A  computer  program  has  been  written  in  BASIC  to  oerform 
the  search  for  the  Lagrange  multiplier  in  the  HMMS  disaggregation 
model.  The  program  has  three  major  parts: 

-  the  code  for  finding  the  bracketing  interval, 

-  the  code  for  oerforming  the  search  with  the  modified 
lattice-Fibonacci  technique, 
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-  the  code  for  Bolzano's  method, 
and  5  subroutines  for: 

-  the  evaluation  of  aggregate  inventory  function  for  a 
given  multiplier  value, 

-  computing  the  derivative  of  the  response  function, 

-  finding  the  asymptote  of  the  response  function, 

-  computing  a  Fibonacci  number  for  the  initial  partitio- 
ning of  the  bracketing  interval, 

-  computing  Fibonacci  numbers  for  a  given  number  of 
experiments. 

The  user  supplies  all  relevant  cost  data,  the  number  of 
items,  the  constraining  aggregate  inventory  and  the  upper  bound 
on  the  error  in  inventory.  The  general  logic  of  the  orogram 
follows  the  development  provided  in  the  paper,  and  is  illustrated 
in  the  attached  flowchart. 

The  program  was  meant  to  serve  the  ouroose  of  the 
comparative  study;  consequently,  the  reports  provide  the  necessary 
information  to  assess  efficiency. 

Some  examples  have  beef!  worked  out  and  tho  reports  are 
shown  below.  A  general  flowchart  as  well  as  the  entire  nrogram 
listing  ere  attached. 


-  47   - 


WHAT    IS    THE    NUMBER    OF     ITEMS? 

!3 

WHAT  15  THE  AGGilEGATE  INVENTORY? 

W}1AT  15  HIGHEST  ACCEPTABLE  ERROR  -  PERCENTAGE  ? 
\P. 


MODIFIED  LATTICEtFIDONACCI  SEARCH  ALGORITHM 

NUMBER  OF  ITEMS  =   3 

IMPOSED  AGGREGATE  INVENTORY  =   800 

MAXIMUM  ADMISSIBLE  ERROR  =■   2  % 

LENGTH  OF  INITIAL  BRACKETING  irOTERVAL  =   .^J41869 

ACTUAL  NUMBER  OF  ITERATIONS  =   5 

ITERATIONS  REQUIRED  3Y  o!JRE  FIBONACCI  SEARCH  =   7 

ESTIMATED  LAGRANGE  MULTIPLIER  =   -.49^^73 

RF.SULTING  AGGREGATE  INVENTORY  =   798.161 

ACTUAL  ERROR  IN  AGGREGATE  INVENTORY  =   -.2298!?8  % 

PRODUCTION  PLAN 

ITEM   1  163  UNITS 

ITEM   2  107  UNITS 

ITEM   3  ^^8  UNITS 

LATTICE  SEARCH  BY  BOLZANO  METHOD 


NUMBER  OF  ITERATIONS  =   4 

ESTIMATED  LAGRANGE  MULTIPLIER  =   -.497103 

RESULTING  AGGREGATE  INVENTORY  =   797.189 

ACTUAL  ERROR  IN  AGGREGATE  INVENTORY  =   -.35136A 

PRODUCTION  PLAN 

ITEM   1         163  UNITS 

ITEM   2        107  UNITS 

ITEM   3        447  UNITS 
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WAT  IS  THE  NUMBER  OF  ITEMS? 

!3 

W^AT  IS  THE  AGGREGATE  INVENTORY? 

WHAT  IS  HIGHEST  ACCEPTABLE  ERROR  -  PERCENTAGE  ? 
!3 


MODIFIED  LATTICE-FIBONACCI  SEARCH  ALGORITHM 


NI'MBER  OF  ITEM<^  =   3 

IM^C-ED  AGGREGA'^E  IN^'ENTOR"  =  ^nci 

MAXIMUM    ADMISSI3LE    ERROR    =       3    % 

LENGTH    OF    INITIAL    BRACKETING    INTER\'AL    =       Rn./(736 

ACTUAL    NUM'lER    OF    ITERATIONS    =        1 '^    , 

ITERATION^-    REQUIRED    3Y    PURE    F,130NACCI     SEARCH    =        13 

ESTIMATED    LAGRANGE    MULTIPLIER    =       -18.1975 

RESULTING    AGGREGATE    INVENTORY    =        199. -963 

ACTUAL    ERROR    IN    AGGREGATE    INVENTORY    =       -1.834  1  lE-HR 

PRODUCTION  PLAN 
ITEM       1  45    UNITS 

ITEM       2  20    UNITS 

ITEM       3  125    UNITS 

LATTICE    SEARCH    BY    BOLZANO    METHOD 


rJUMBER    OFMTERATIONS    =       9 

ESTIMATED    LAGRANGE    MULTIPLIER    =       -18.1966 

RESULTING    AGGREGATE    INVENTORY    =        199.968 

ACTUAL    ERROR    IN    AGGREGATE    INVENTORY    =       -  1  .  6 1  ?  85E-'?!2    % 

PRODUCTION  PLAN 
ITEM       1  4  5    UNITS 

ITEM       -5  20    UNITS 

ITEM      3  125    UNITS 
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APPENDIX  I 


Subroutine  3: 
find  asymptote 


^„  =  min 


'11 


Subroutine  2; 
comoute 


din 


Subroutine  1; 
compute  Iq( A ) 


Input  number  of 
items,  aggregate 
inventory,  and 
bound  on  the  error 


Bracket  the  inter- 
val containing  the 
Lagrange  multiolier 


Partition  the 
bracketing  inter- 
val 


Perform  search 
for  multiplier 


Print  report  from 
lattice-Fibonacci 
search 


Perform  Bolzano 
search  for  multipl. 


Print  report  from 
Bolzano  search 


Subroutine  4*. 
compute  Fibonacci 
number  for 
initial  parti- 
tioning 


Subroutine  5: 
compute  Fibonacci 
number  for 
given  number 
of  experiments 


■Flow  chart  ahowing  the  general  features  of  the 
cQiroarative  search  program 
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APPENDIX  II  :  Comouter  orogram  for  the  modified  lattice-Fibonacci 
search  and  for  Bolzano's  root  finding  technique 


:10 

20 

30 

40 

50 

f.Pi 

70 

80 

90 

in0 

1:1  0 

±20 

l-:0 

140 

:1.50 

160 

170 

180 

±90 

200 

210 
220 
220 
240 
250 
260 
270 
280 
290 
j:00 
310 
320 
330 
340 
3!r.0 
360 
370 
380 
390 
480 
4J0 
420 
4~:0 
440 
450 
4i:.H 
470 
480 
490 
300 
310 
320 
330 


PRINI   WHfiT  15  THE  NUMBER  OF  ITEMS?- 

INPUT  I 

FOR  K>=i  TO  I 

READ  c<K>..  h<::k:>..  U':;K>..  S<hO 

NEXT  K 

PRINI  'WHRT  IS  THE  RQGREGRTE  INVENTORV?' 

INPUT  R 

LET  W2--^l 

LET  X^0 

GOSUB  1430 

IF  RG-=R  THEN  1230 

GOSUB  1490 

REM  +  FIND  INTERSECTION  OF  TRNGENT  RND  HORIZONTAL  RXIS 

LET  K-^<R-R0>.-'D 

REM  *  FIND  BRRCKETING  INTERVRL  * 

IF  X<0  THEN  340 

GOSUB  1330 

IF  K::=-M  THEN  230 

LET  E:1^0 

LET  E2=X 

GOTO  600 

REI'1  +  BINRRV  SERRCH  FOR  FINDING  BRRCKETING  INTERVRL  + 

LET  E1^0 

LET  E2=^El+'::M-Ei::'/2 

LET  X=E2 

GOSUB  1430 

LET  V^R0-R 

IF  V-G  THEN  1230 

IF  V<:0  THEN  310 

GOTO  430 

LET  E1--^E2 

GOTO  240 

RFM  *  BINRRV  SERRCH  ENDS  HERE  + 

LET  E2=X 

LET  K>1 

LET  L-^2-K*E2 

LET  X"-=,X+L 

GOSUB  1430 

LET  V^R0-R 

IF  V-^0  THEN  1270 

IF  VCO  THEN  440 

LET  K^K+1 

GOTO  360 

LET  E1=X 

REM  +■+  PRRTITION  BRRCKETIflG  INTERVRL  -::E1.  E2>  ** 

PRINT  'WHRT  IS  HIGHEST  RCC EPTRBLE  ERROR  -  PERCENTRGE  ? 

INPUT  B 

LET  X^E2 

GOSUB  1430 

GOSUB  1490 

LET  T^^D* ''  E2-E1  :■  /  <  RG  +  B/lOO 

LET  E3-^E1 

LET  E4-E2 
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540 
"550 
56  0 
570 
590 
5f-'0 

600 

610 
6?0 
670 
6'10 

6i:-0 
€7  id 
6  90 
6l-'0 

700 

710 
720 
7~'0 
740 
750 
760 
770 
780 
790 

900 

8-1  0 
820 
870 
840 
8r'0 
8^-0 
870 
880 
890 

900 

910 
920 
930 
940 
950 
960 
970 

990 

1000 
10J  0 
1020 
1070 
1040 
1050 
1060 


810 


GOSUB  2260 

LET  L-E2-E1 

GOSI.IB  1660 

LET  N0-=N 

LET  W=^0 

Fi-EM  +■++  MODIFIED  LRTTICE-FIBONRCCI 

LET  Xi^-El 

LET  X--^X1 

GOSUB  1430 

LET  V1-FI0-R 

LET  X2^-^E? 

LET  Jv^X2 

GOSUB  1430 

LET  V2-R0-R 

IF  F-2  THEN  ."10 

IF  F>2  THEN  830 

GOTO  1130 

LET  X3~<Xi+X2:>/2 

LET  X=.X3 

GOSUB  1430 

LET  W^W+1 

V3~R0-R 

IF  V3=0  THEN  1290 

LET  F-^1 

IF  V3::0  THEN 

LET  E2==X3 

GOTO  600 

LET  E1=X3 

GOTO  600 

LET  .T--^N-1 

GOSUB  1780 

LET  X3=X1+G+ 

LET  y.^K2 

GOSUB  1430 

LET  N-^W+1 

LET  V3=R0-R 

IF  V3^0  THEN  129f 

LET  X4^^:-^;2-G+ 

LET  :v^X4 

GOSUB  1430 

LET  N--^W+1 

LFT  V4=R0-R 

IF  V4=^0  THEN  1310 

IF  V4<0  THEN  SS<d 

GOTO  1010 

IF  V3<0  THEN  1070 

GOTO    lHi9iZ-i 

LET  E2-X4 

LET  N=N-2 

LET  T-N 

GOSUB  1780 

LET  F-G 

GOTO  600 


^.ERRCH  RLGORITHM  ++•* 


E2-Ei:>.'F 


E2-Ei:>.'F 
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LET  E1^X3 

GOTO  J.028 

LFT  Ei-X4 

LET  E2=^X3 

LET  r^=N-3 

GOTO  1030 

REM  ++THE  INTERVRL  IS  1  UNIT  LONG  RND  CONTRINS  NO  POINT++ 

REM  *■■    INTERVfll.  END  WITH  THE  SMRLLEST  ERROR  IN  + 

REM  +■  HGGREGflTE  INVENTORV  IS  MULTIRLIER  ESTIMATE  +• 

LET  Vl^RBS-iiVl;- 

IF  V1<:V2  THEN  1210 

REM  +  3ERRCH  IS  COMPLETE.;  X0  IS  MULTIPLIER  ESTIMRTE  + 


LET  ? 

S0-X2 

GOTO 

1320 

LET  J 

^0^X1 

GOTO 

132© 

LET  :■ 

•••:0-0 

GOTO 

1320 

LET  ) 

^0=E2 

GOTO 

1320 

LET  ; 

•^0-:k 

■'R  ■•+100 


:k>*Svk 


>/<H':;k:>-x+U'::k> 


GOTO  1320 

LET  X0=K3 

GOTO  1320 

LET  X0-K4 

LET  X-X0 

GOSUB  1430 

LET  P-''. '::R0-R::' 

FOR  K^l  TO  I 

LEI  Q<K>^'::2+i: 

LET  Q-:;K>=^INT':; 

NEXT  K 

IF  W2-2  THFN  2050 

GOTO  1920 

REM    ++■••+•+:+    SUBROUTINES 

REM  ++  SUBROUTINE  FOR 

LET  R0-0 

FOR  K-=l  TO  I 

LET  RQ-fie+  •:; 2+C •:; K  >  +S  •::  K 

rJEXT  K 

RETURN 

REM  +•+  SUBROUTINE  FOR  COMfUTING  FIRST  DERIVRTIVE  *+ 

LET  D--^0 

FOR  K>1  TO  I 

LEI  D^D+'::2+C^K::'+S'::K 

NEXT  K 

RETURN 

PFMi+SUBROUTINE  FOR  FINDItiG  MINIMUM  H/U** 


COMPUTING 


RGQREQRTE  INVENTORV  +* 


:h-  K 


-X*U  ■:;  K  ::•>>■".  5*U  <  K  >  .•■'2 


;4*<H':;k:;'-x*U''K>>-i. 


FOR    K:-=l 

LET  R'.:k; 

r-JEXT    K 
LET    M^-R> 
FOR    K^2 


TO    I 
^H'::K 


1' 
TO 


■ij  v  K : 
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IF  PvK>  :M  then  1620 

GOTO  1620 

LET  M-R'::K> 

NEXT  t 

RETUFTJ 

F:En  +-+ SUBROUTINE  1  FOR 

FIBi 

IiNflC 

C I  NUMBERS** 

LET  F0-i 

LET  Fl-^i 

LET  N-1 

LET  F-Fl 

IF  F::--T  then  17613 

LET  F--Fi+F0 

LET  F0^:F1 

LET  Fl-F 

LET  N-^N+1 

GOTO  1700 

RETURrJ 

RFI'1+  + SUBROUTINE  2  FOR 

FIBONHCC 

:  I  NUMBERS** 

IF  J-0  THEN  lE.'iO 

IF  J-1  THEN  1S10 

GOTO  1830 

LET  G-1 

GOTO  1900 

LET  G0--^i 

LET  Gl-i 

FOR  K-^2  TO  J 

LET  G-G1+G0 

LET  G0--^G1 

LET  Gl^G 

NE:a  K 

HERE***** 


MODIFIEC:-    LFlTTICE-FIBONflCCI    SEARCH    RLGORITHM' 


RETURN 

REM  *****RLL  SUBROUTINES  END 

DflTfl  10..  1..  1..  2000..  10..  1..  5..  2000 

DBTfi  30..  2.  2..  10000 

PRINT 

PRINT 

PRINT 

PRINT 

PR  I  NT 

PRINT 

PRINT 

PRINT 

PRINT 

PRINT 

IF  l.J2=l  THEN  2030 

PRINT 

PRTNT  ••LflTTICE  SEARCH  BV  BOLZANO  METHOD' 

PPT  NT   '^^^.^::^^:^^:^^:^:^^::^^:=^==^^=--^  =  =  ^^:^^^==^-- 

PPINT 

IF  W2-2  THEN  2140 

PRINT  'LENGTH  OF  INITIAL  E^:ACKETING  INTERVAL  =  '  :  L 

PRINT   ACTUAL  NUMBER  OF  HIRATIONS  ■■=     ':W 

PRINT  'ITERATIONS  REQUIREC  BV  PURE  FIBONACCI  SEARCH 


NUMBER  OF  ITEMS  ^  ■  ; I 

IMPOSED  AGGREGATE  INVENTORV  =  • 

MAXIMUM  ADMISSIBLE  ERROR  ■-  ';B 


A 


N0 
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IF  W2-^l  THEN  2150 

PRINT  'NUMBER  OF  ITERATIONS  -=   :  Wi 

PRINT 

PRINT  •ESTIMFITED  LFlQRRNQE  MULTIPLIER  ■=   :X 

PRINT  'RESULTING  RQGREGFlTE  INVENTORY  -^    •-  :  09 

PRINT  -FlCTUfiL  ERROR  IN  RGGREGFlTE  INVENTORY  =  ';  P  ;:•••; •- 

PRINT 

PRINT  'PRODUCTION  PLRN' 

FOR  K-^1  TO  I 

PR  TNT  'ITEM  '  :  K .  Q  •'  K  >  :  -  UN  I TS  •• 

NEXT  K 

GOTO  21' 30 

RFM  +*  PRRTITION  BRRCKETING  INTERVRL  FOR  BOLZRNO  SERRCH  *+ 

LFT  K>1 

LET  V^-=2-K 

IF  V>=^T  THEN  2310 

LET  K^K+l 

GOTO  2276 

RETURN 

REM  +*  BOLZRNO  SERRCH  ++• 

IF  W2^2  THEN  2639 

LET  wi-=e 

LET  Vl-^V/2 

LET  X^=^E3+V1+<E4-E3;>/V 

GOSUB  1420 

LET  W1=W1+1 

LFT  V-^RO-R 

LFT  XO=-X 

IF  V-^-e  THEN  13:20 

ir  V<0  THEN  2450 

LFT  E4-=X 

GOTO  2460 

LET  E3:=X 

IF  V1--1  THEN  2490 

LET  V-V-Vi 

GOTO  2350 

LET  v=E3: 

GOSUB  1420 

LET  V1-^R0-R 

LFT  Yl^RBS'::Vi::' 

LET  X=E4 

GOSUB  1420 

LFT  V2=^R0-R 

IF    VK-V2    THEN    2590 

LET    X0-^E4 

GriTO   2600 

LET    X0--==E3 

LFT    N2=^2 

GOTO    1320 

RFM    +•+••    END    OF    BOLZRNO    SERF  CH    *+ 

EfJD 


